#################################################
### Replication Code for Study 1              ###
### Title: Mimicking Democracy                ###
### Authors: Hsu Yumin Wang, Eddy S. F. Yeung ###
### Version: April 26, 2024                   ###
#################################################

### Set-up ----
rm(list = ls())
setwd("~/Desktop/redist-propaganda/replication/study 1")
library(tidyverse)   # version 2.0.0
library(survey)      # version 1.1-1
library(cowplot)     # version 1.1.1
library(ggrepel)     # version 0.9.1
library(grid)        # version 4.0.1
library(gridExtra)   # version 2.3
library(lme4)        # version 1.1-30
library(lmerTest)    # version 3.1-3
library(arm)         # version 1.13-1
library(estimatr)    # version 1.0.0
library(texreg)      # version 1.37.5
library(countrycode) # version 1.3.0
library(haven)       # version 2.4.3
library(mice)        # version 3.14.0

### Import datasets ----
df <- read.csv("GBS3_20200616.csv")
country_label <- read.csv("country_label_study1.csv")
nrow(df)
df <- subset(df, wt > 0) # drop if weight is negative
nrow(df)

### Recode variables ----
## DV: democratic satisfaction (0 = not at all satisfied; 3 = very satisfied)
table(df$sd1)
df <- df %>% 
  mutate(dem_satis = case_when(
    sd1 == 0            ~ 0,
    sd1 >= 1 & sd1 <= 4 ~ sd1 - 1
  ))
table(df$sd1, df$dem_satis)

## DV: democratic evaluation (0 = not a democracy; 3 = a full democracy)
table(df$sd2)
df <- df %>% 
  mutate(dem_eval = case_when(
    sd2 >= 1 & sd2 <= 4 ~ sd2 - 1
  ))
table(df$sd2, df$dem_eval)

## IV: govt performance in narrowing income gaps (0 = very bad; 3 = very good)
table(df$pcg3)
df <- df %>% 
  mutate(perf_incgap = case_when(
    pcg3 >= 1 & pcg3 <= 4 ~ pcg3 - 1
  ))
table(df$pcg3, df$perf_incgap)

## IV: perceived fairness of income distribution (0 = very unfair; 3 = very fair)
table(df$rd1)
df <- df %>% 
  mutate(fair_incdist = case_when(
    rd1 >= 1 & rd1 <= 4 ~ rd1 - 1
  ))
table(df$rd1, df$fair_incdist)

## Control: govt performance in managing the economy (0 = very bad; 3 = very good)
table(df$pcg1)
df <- df %>% 
  mutate(perf_econ = case_when(
    pcg1 >= 1 & pcg1 <= 4 ~ pcg1 - 1
  ))
table(df$pcg1, df$perf_econ)

## Control: govt performance in creating jobs (0 = very bad; 3 = very good)
table(df$pcg2)
df <- df %>% 
  mutate(perf_job = case_when(
    pcg2 >= 1 & pcg2 <= 4 ~ pcg2 - 1
  ))
table(df$pcg2, df$perf_job)

## Control: govt performance in keeping prices down (0 = very bad; 3 = very good)
table(df$pcg4)
df <- df %>% 
  mutate(perf_price = case_when(
    pcg4 >= 1 & pcg4 <= 4 ~ pcg4 - 1
  ))
table(df$pcg4, df$perf_price)

## Control: govt performance in providing health services (0 = very bad; 3 = very good)
table(df$pcg5)
df <- df %>% 
  mutate(perf_health = case_when(
    pcg5 >= 1 & pcg5 <= 4 ~ pcg5 - 1
  ))
table(df$pcg5, df$perf_health)

## Control: govt performance in addressing educational needs (0 = very bad; 3 = very good)
table(df$pcg6)
df <- df %>% 
  mutate(perf_educ = case_when(
    pcg6 >= 1 & pcg6 <= 4 ~ pcg6 - 1
  ))
table(df$pcg6, df$perf_educ)

## Demographics: gender (0 = male; 1 = female)
table(df$gbse1)
df <- df %>% 
  mutate(female = case_when(
    gbse1 == 1 ~ 0,
    gbse1 == 2 ~ 1
  ))
table(df$gbse1, df$female)

## Demographics: age (16 = minimum; 108 = maximum)
table(df$gbse2)
df <- df %>% 
  mutate(age = case_when(
    gbse2 != -1 ~ gbse2
  ))

## Demographics: education (0 = primary or no education; 4 = university or above)
table(df$gbse3)
df <- df %>% 
  mutate(educ = case_when(
    gbse3 >= 1 & gbse3 <= 5 ~ gbse3 - 1
  ))
table(df$gbse3, df$educ)

## Demographics: Internet usage (0 = never; 3 = daily)
table(df$int1)
df <- df %>% 
  mutate(internet = case_when(
    int1 == 0             ~ 0,
    int1 >= 1 & int1 <= 4 ~ int1 - 1
  ))
table(df$int1, df$internet)

### Cross-country analysis ----
## Summarize mean DVs and IVs by country by applying survey weights
df_weight <- svydesign(ids = ~0, strata = NULL, weights = ~wt, data = df)
df_dem_satis <- svyby(~dem_satis, ~country_un, df_weight, svymean, na.rm = T)
df_dem_eval <- svyby(~dem_eval, ~country_un, df_weight, svymean, na.rm = T)
df_perf_incgap <- svyby(~perf_incgap, ~country_un, df_weight, svymean, na.rm = T)
df_fair_incdist <- svyby(~fair_incdist, ~country_un, df_weight, svymean, na.rm = T)
df_summary <- country_label %>% 
  left_join(df_dem_satis, by = "country_un") %>% 
  left_join(df_dem_eval, by = "country_un") %>% 
  left_join(df_perf_incgap, by = "country_un") %>% 
  left_join(df_fair_incdist, by = "country_un")
df_summary$dem_satis <- ifelse(df_summary$dem_satis == 0, NA, df_summary$dem_satis)
df_summary$dem_eval <- ifelse(df_summary$dem_eval == 0, NA, df_summary$dem_eval)
df_summary$perf_incgap <- ifelse(df_summary$perf_incgap == 0, NA, df_summary$perf_incgap)
df_summary$fair_incdist <- ifelse(df_summary$fair_incdist == 0, NA, df_summary$fair_incdist)

## Democratic satisfaction and govt performance in narrowing income gaps
p1 <- 
  ggplot(df_summary, aes(x = perf_incgap, y = dem_satis), label = country_label) + 
  geom_point(shape = 20, size = 2) + 
  stat_smooth(method = "loess", formula = y ~ x, se = F, na.rm = T,  
              fullrange = T, linetype = 2, color = "blue") + 
  geom_text_repel(aes(x = perf_incgap, y = dem_satis, label = country_label), 
                  size = 2.5, max.overlaps = 12, family = "Times") + 
  coord_cartesian(xlim = c(0, 2.5), ylim = c(0, 2.5)) + 
  labs(x = "Government Performance in Narrowing Income Gaps\n(0 = very bad; 3 = very good)", 
       y = "Democratic Satisfaction\n(0 = not at all satisfied; 3 = very satisfied)") + 
  theme_bw() + 
  ggtitle("Afrobarometer (Wave 6) & Arab Barometer (Wave 4)") + 
  theme(text = element_text(family = "Times", size = 12), 
        axis.text = element_text(color = "black"), 
        plot.title = element_text(size = 12, face = "italic"))
r <- cor(x = df_summary$perf_incgap, y = df_summary$dem_satis, use = "complete.obs") # correlation coefficient
n <- sum(!is.na(df_summary$perf_incgap) & !is.na(df_summary$dem_satis)) # count the number of complete cases
SE <- sqrt((1 - r^2) / (n - 2)) # SE of the correlation coefficient
p1 <- p1 +
  annotate("label", x = 2.25, y = 2.25, family = "Times",
           label = paste0("r = ", format(r, digits = 2), "\n(SE = ", format(SE, digits = 2), ")"), size = 2.5)

## Democratic evaluation and govt performance on narrowing income gaps
p2 <- 
  ggplot(df_summary, aes(x = perf_incgap, y = dem_eval), label = country_label) + 
  geom_point(shape = 20, size = 2) + 
  stat_smooth(method = "loess", formula = y ~ x, se = F, na.rm = T,  
              fullrange = T, linetype = 2, color = "blue") + 
  geom_text_repel(aes(x = perf_incgap, y = dem_eval, label = country_label), 
                  size = 2.5, max.overlaps = 12, family = "Times") + 
  coord_cartesian(xlim = c(0, 2.5), ylim = c(0, 2.5)) + 
  labs(x = "Government Performance in Narrowing Income Gaps\n(0 = very bad; 3 = very good)", 
       y = "Democratic Evaluation\n(0 = not a democracy; 3 = a full democracy)") + 
  theme_bw() + 
  ggtitle("Afrobarometer (Wave 6) & Arab Barometer (Wave 4)") + 
  theme(text = element_text(family = "Times", size = 12), 
        axis.text = element_text(color = "black"), 
        plot.title = element_text(size = 12, face = "italic"))
r <- cor(x = df_summary$perf_incgap, y = df_summary$dem_eval, use = "complete.obs") # correlation coefficient
n <- sum(!is.na(df_summary$perf_incgap) & !is.na(df_summary$dem_eval)) # count the number of complete cases
SE <- sqrt((1 - r^2) / (n - 2)) # SE of the correlation coefficient
p2 <- p2 +
  annotate("label", x = 2.25, y = 2.25, family = "Times",
           label = paste0("r = ", format(r, digits = 2), "\n(SE = ", format(SE, digits = 2), ")"), size = 2.5)

## Democratic satisfaction and perceived fairness of income distribution
p3 <- 
  ggplot(df_summary, aes(x = fair_incdist, y = dem_satis), label = country_label) + 
  geom_point(shape = 20, size = 2) + 
  stat_smooth(method = "loess", formula = y ~ x, se = F, na.rm = T,  
              fullrange = T, linetype = 2, color = "blue") + 
  geom_text_repel(aes(x = fair_incdist, y = dem_satis, label = country_label), 
                  size = 2.5, max.overlaps = 12, family = "Times") + 
  coord_cartesian(xlim = c(0, 2.5), ylim = c(0, 2.5)) + 
  labs(x = "Perceived Fairness of Income Distribution\n(0 = very unfair; 3 = very fair)", 
       y = "Democratic Satisfaction\n(0 = not at all satisfied; 3 = very satisfied)") + 
  theme_bw() + 
  ggtitle("Asian Barometer (Wave 4) & Latinobarómetro (Wave 15)") + 
  theme(text = element_text(family = "Times", size = 12), 
        axis.text = element_text(color = "black"), 
        plot.title = element_text(size = 12, face = "italic"))
r <- cor(x = df_summary$fair_incdist, y = df_summary$dem_satis, use = "complete.obs") # correlation coefficient
n <- sum(!is.na(df_summary$fair_incdist) & !is.na(df_summary$dem_satis)) # count the number of complete cases
SE <- sqrt((1 - r^2) / (n - 2)) # SE of the correlation coefficient
p3 <- p3 +
  annotate("label", x = 2.25, y = 2.25, family = "Times",
           label = paste0("r = ", format(r, digits = 2), "\n(SE = ", format(SE, digits = 2), ")"), size = 2.5)

## Democratic evaluation and perceived fairness of income distribution
p4 <- 
  ggplot(df_summary, aes(x = fair_incdist, y = dem_eval), label = country_label) + 
  geom_point(shape = 20, size = 2) + 
  stat_smooth(method = "loess", formula = y ~ x, se = F, na.rm = T,  
              fullrange = T, linetype = 2, color = "blue") + 
  geom_text_repel(aes(x = fair_incdist, y = dem_eval, label = country_label), 
                  size = 2.5, max.overlaps = 12, family = "Times") + 
  coord_cartesian(xlim = c(0, 2.5), ylim = c(0, 2.5)) + 
  labs(x = "Perceived Fairness of Income Distribution\n(0 = very unfair; 3 = very fair)", 
       y = "Democratic Evaluation\n(0 = not a democracy; 3 = a full democracy)") + 
  theme_bw() + 
  ggtitle("Asian Barometer (Wave 4)") + 
  theme(text = element_text(family = "Times", size = 12), 
        axis.text = element_text(color = "black"), 
        plot.title = element_text(size = 12, face = "italic"))
r <- cor(x = df_summary$fair_incdist, y = df_summary$dem_eval, use = "complete.obs") # correlation coefficient
n <- sum(!is.na(df_summary$fair_incdist) & !is.na(df_summary$dem_eval)) # count the number of complete cases
SE <- sqrt((1 - r^2) / (n - 2)) # SE of the correlation coefficient
p4 <- p4 +
  annotate("label", x = 2.25, y = 2.25, family = "Times",
           label = paste0("r = ", format(r, digits = 2), "\n(SE = ", format(SE, digits = 2), ")"), size = 2.5)

# Combine into one graph (Figure S3)
cross_country_cor <- plot_grid(p1, p2, p3, p4, labels = "AUTO", ncol = 2, 
                               label_fontfamily = "Times")
ggsave(file = "cross_country_cor.pdf", cross_country_cor, width = 9, height = 9)

## Placebo tests based on actual levels of democracy
# Merge with V-Dem's Electoral Democracy index
vdem <- read_dta("vdem.dta")
vdem$country_un <- countrycode(vdem$cowcode, origin = 'cown', destination = 'un') # create country_un in V-Dem for merging
vdem$year_lag <- vdem$year + 1 # V-Dem index (t-1) and survey data (t)
df$year_lag <- df$year
df <- left_join(df, vdem, by = c("country_un", "year_lag"))
df_vdem <- df %>% 
  group_by(country_un) %>% 
  summarize(mean_v2x_polyarchy = mean(v2x_polyarchy))
df_summary <- df_summary %>% 
  left_join(df_vdem, by = "country_un")

# Electoral Democracy index and govt performance in narrowing income gaps
p5 <- 
  ggplot(df_summary, aes(x = perf_incgap, y = mean_v2x_polyarchy), label = country_label) + 
  geom_point(shape = 20, size = 2) + 
  stat_smooth(method = "loess", formula = y ~ x, se = F, na.rm = T,  
              fullrange = T, linetype = 2, color = "blue") + 
  geom_text_repel(aes(x = perf_incgap, y = mean_v2x_polyarchy, label = country_label), 
                  size = 2.5, max.overlaps = 12, family = "Times") + 
  coord_cartesian(xlim = c(0, 2.5), ylim = c(0, 1)) + 
  labs(x = "Government Performance in Narrowing Income Gaps\n(0 = very bad; 3 = very good)", 
       y = "V-Dem Index (v2x_polyarchy)") + 
  theme_bw() + 
  ggtitle("Afrobarometer (Wave 6) & Arab Barometer (Wave 4)") + 
  theme(text = element_text(family = "Times", size = 12), 
        axis.text = element_text(color = "black"), 
        plot.title = element_text(size = 12, face = "italic"))
r <- cor(x = df_summary$perf_incgap, y = df_summary$mean_v2x_polyarchy, use = "complete.obs") # correlation coefficient
n <- sum(!is.na(df_summary$perf_incgap) & !is.na(df_summary$mean_v2x_polyarchy)) # count the number of complete cases
SE <- sqrt((1 - r^2) / (n - 2)) # SE of the correlation coefficient
p5 <- p5 +
  annotate("label", x = 2.25, y = .875, family = "Times",
           label = paste0("r = ", format(r, digits = 2), "\n(SE = ", format(SE, digits = 2), ")"), size = 2.5)

# Electoral Democracy index and perceived fairness of income distribution
p6 <- 
  ggplot(df_summary, aes(x = fair_incdist, y = mean_v2x_polyarchy), label = country_label) + 
  geom_point(shape = 20, size = 2) + 
  stat_smooth(method = "loess", formula = y ~ x, se = F, na.rm = T,  
              fullrange = T, linetype = 2, color = "blue") + 
  geom_text_repel(aes(x = fair_incdist, y = mean_v2x_polyarchy, label = country_label), 
                  size = 2.5, max.overlaps = 12, family = "Times") + 
  coord_cartesian(xlim = c(0, 2.5), ylim = c(0, 1)) + 
  labs(x = "Perceived Fairness of Income Distribution\n(0 = very unfair; 3 = very fair)", 
       y = "V-Dem Index (v2x_polyarchy)") + 
  theme_bw() + 
  ggtitle("Asian Barometer (Wave 4) & Latinobarómetro (Wave 15)") + 
  theme(text = element_text(family = "Times", size = 12), 
        axis.text = element_text(color = "black"), 
        plot.title = element_text(size = 12, face = "italic"))
r <- cor(x = df_summary$fair_incdist, y = df_summary$mean_v2x_polyarchy, use = "complete.obs") # correlation coefficient
n <- sum(!is.na(df_summary$fair_incdist) & !is.na(df_summary$mean_v2x_polyarchy)) # count the number of complete cases
SE <- sqrt((1 - r^2) / (n - 2)) # SE of the correlation coefficient
p6 <- p6 +
  annotate("label", x = 2.25, y = .875, family = "Times",
           label = paste0("r = ", format(r, digits = 2), "\n(SE = ", format(SE, digits = 2), ")"), size = 2.5)

# Combine into one graph (Figure S4)
cross_country_placebo <- plot_grid(p5, p6, labels = "AUTO", ncol = 2, 
                                   label_fontfamily = "Times")
ggsave(file = "cross_country_placebo.pdf", cross_country_placebo, width = 9, height = 4.5)

### Within-country analysis (multilevel modeling) ----
## Democratic satisfaction and govt performance in narrowing income gaps
# Impute the missing data
df$country_un <- as.factor(df$country_un)
table(df$pcg3, df$country_un) # identify countries that were asked the IV question
df2 <- df %>%
  dplyr::filter(country_un == 12 | country_un == 72 | country_un == 108 |
                  country_un == 120 | country_un == 132 | country_un == 204 |
                  country_un == 266 | country_un == 275 | country_un == 288 |
                  country_un == 324 | country_un == 384 | country_un == 400 |
                  country_un == 404 | country_un == 422 | country_un == 426 |
                  country_un == 430 | country_un == 450 | country_un == 454 |
                  country_un == 466 | country_un == 480 | country_un == 504 |
                  country_un == 508 | country_un == 516 | country_un == 562 |
                  country_un == 566 | country_un == 678 | country_un == 686 |
                  country_un == 694 | country_un == 710 | country_un == 716 |
                  country_un == 748 | country_un == 768 | country_un == 788 |
                  country_un == 800 | country_un == 818 | country_un == 834 |
                  country_un == 854 | country_un == 894) %>%
  dplyr::select(country_un, dem_satis, perf_incgap,
                female, age, educ, internet, perf_econ, perf_job,
                perf_price, perf_health, perf_educ) %>% 
  dplyr::filter(!is.na(dem_satis)) # drop respondents with missing DV
imputed_data <- mice(df2, seed = 1234567)
df_hlm_A <- complete(imputed_data, 1)
mean(df_hlm_A$dem_satis) # baseline mean of the DV = 1.42
sd(df_hlm_A$dem_satis) # standard deviation of the DV = 0.99

# Estimate the multilevel models
hlm_A1 <- lmer(dem_satis ~ perf_incgap + 
                 (1 + perf_incgap | country_un), 
               data = df_hlm_A,
               control = lmerControl(optimizer = "Nelder_Mead"))
hlm_A2 <- lmer(dem_satis ~ perf_incgap + 
                 female + age + I(age^2/100) + educ + internet + 
                 (1 + perf_incgap | country_un), 
               data = df_hlm_A,
               control = lmerControl(optimizer = "Nelder_Mead"))
hlm_A3 <- lmer(dem_satis ~ perf_incgap + 
                 female + age + I(age^2/100) + educ + internet + 
                 perf_econ + perf_job + perf_price + perf_health + perf_educ + 
                 (1 + perf_incgap | country_un), 
               data = df_hlm_A,
               control = lmerControl(optimizer = "Nelder_Mead"))

## Democratic evaluation and govt performance in narrowing income gaps
# Impute the missing data
df3 <- df %>%
  dplyr::filter(country_un == 12 | country_un == 72 | country_un == 108 |
                  country_un == 120 | country_un == 132 | country_un == 204 |
                  country_un == 266 | country_un == 275 | country_un == 288 |
                  country_un == 324 | country_un == 384 | country_un == 400 |
                  country_un == 404 | country_un == 422 | country_un == 426 |
                  country_un == 430 | country_un == 450 | country_un == 454 |
                  country_un == 466 | country_un == 480 | country_un == 504 |
                  country_un == 508 | country_un == 516 | country_un == 562 |
                  country_un == 566 | country_un == 678 | country_un == 686 |
                  country_un == 694 | country_un == 710 | country_un == 716 |
                  country_un == 748 | country_un == 768 | country_un == 788 |
                  country_un == 800 | country_un == 818 | country_un == 834 |
                  country_un == 854 | country_un == 894) %>%
  dplyr::select(country_un, dem_eval, perf_incgap,
                female, age, educ, internet, perf_econ, perf_job,
                perf_price, perf_health, perf_educ) %>% 
  dplyr::filter(!is.na(dem_eval)) # drop respondents with missing DV
imputed_data <- mice(df3, seed = 1234567)
df_hlm_B <- complete(imputed_data, 1)
mean(df_hlm_B$dem_eval) # baseline mean of the DV = 1.62
sd(df_hlm_B$dem_eval) # standard deviation of the DV = 0.91

# Estimate the multilevel models
hlm_B1 <- lmer(dem_eval ~ perf_incgap + 
                 (1 + perf_incgap | country_un), 
               data = df_hlm_B,
               control = lmerControl(optimizer = "Nelder_Mead"))
hlm_B2 <- lmer(dem_eval ~ perf_incgap + 
                 female + age + I(age^2/100) + educ + internet + 
                 (1 + perf_incgap | country_un), 
               data = df_hlm_B,
               control = lmerControl(optimizer = "Nelder_Mead"))
hlm_B3 <- lmer(dem_eval ~ perf_incgap + 
                 female + age + I(age^2/100) + educ + internet + 
                 perf_econ + perf_job + perf_price + perf_health + perf_educ + 
                 (1 + perf_incgap | country_un), 
               data = df_hlm_B, 
               control = lmerControl(optimizer = "Nelder_Mead"))

## Democratic satisfaction and perceived fairness of income distribution
# Impute the missing data
table(df$rd1, df$country_un) # identify countries that were asked the IV question
df4 <- df %>%
  dplyr::filter(country_un == 32 | country_un == 68 | country_un == 76 |
                  country_un == 104 | country_un == 116 | country_un == 152 |
                  country_un == 156 | country_un == 170 | country_un == 188 |
                  country_un == 214 | country_un == 218 | country_un == 222 |
                  country_un == 320 | country_un == 340 | country_un == 344 |
                  country_un == 360 | country_un == 392 | country_un == 410 |
                  country_un == 458 | country_un == 484 | country_un == 496 |
                  country_un == 558 | country_un == 591 | country_un == 600 |
                  country_un == 604 | country_un == 608 | country_un == 702 |
                  country_un == 704 | country_un == 713 | country_un == 764 |
                  country_un == 858 | country_un == 862) %>%
  dplyr::select(country_un, dem_satis, fair_incdist,
                female, age, educ, internet) %>% 
  dplyr::filter(!is.na(dem_satis)) # drop respondents with missing DV
imputed_data <- mice(df4, seed = 1234567)
df_hlm_C <- complete(imputed_data, 1)
mean(df_hlm_C$dem_satis) # baseline mean of the DV = 1.53
sd(df_hlm_C$dem_satis) # standard deviation of the DV = 0.84

# Estimate the multilevel models
hlm_C1 <- lmer(dem_satis ~ fair_incdist + 
                 (1 + fair_incdist | country_un),
               data = df_hlm_C,
               control = lmerControl(optimizer = "Nelder_Mead"))
hlm_C2 <- lmer(dem_satis ~ fair_incdist + 
                 female + age + I(age^2/100) + educ + internet + 
                 (1 + fair_incdist | country_un), 
               data = df_hlm_C,
               control = lmerControl(optimizer = "Nelder_Mead"))

## Democratic evaluation and perceived fairness of income distribution
# Impute the missing data
df5 <- df %>%
  dplyr::filter(country_un == 32 | country_un == 68 | country_un == 76 |
                  country_un == 104 | country_un == 116 | country_un == 152 |
                  country_un == 156 | country_un == 170 | country_un == 188 |
                  country_un == 214 | country_un == 218 | country_un == 222 |
                  country_un == 320 | country_un == 340 | country_un == 344 |
                  country_un == 360 | country_un == 392 | country_un == 410 |
                  country_un == 458 | country_un == 484 | country_un == 496 |
                  country_un == 558 | country_un == 591 | country_un == 600 |
                  country_un == 604 | country_un == 608 | country_un == 702 |
                  country_un == 704 | country_un == 713 | country_un == 764 |
                  country_un == 858 | country_un == 862) %>%
  dplyr::select(country_un, dem_eval, fair_incdist,
                female, age, educ, internet) %>% 
  dplyr::filter(!is.na(dem_eval)) # drop respondents with missing DV
imputed_data <- mice(df5, seed = 1234567)
df_hlm_D <- complete(imputed_data, 1)
mean(df_hlm_D$dem_eval) # baseline mean of the DV = 1.73
sd(df_hlm_D$dem_eval) # standard deviation of the DV = 0.75

# Estimate the multilevel models
hlm_D1 <- lmer(dem_eval ~ fair_incdist + 
                 (1 + fair_incdist | country_un), 
               data = df_hlm_D,
               control = lmerControl(optimizer = "Nelder_Mead"))
hlm_D2 <- lmer(dem_eval ~ fair_incdist + 
                 female + age + I(age^2/100) + educ + internet + 
                 (1 + fair_incdist | country_un), 
               data = df_hlm_D,
               control = lmerControl(optimizer = "Nelder_Mead"))

## Table 1: Panels A and B (mixed effects)
texreg(list(hlm_A1, hlm_A2, hlm_A3, hlm_B1, hlm_B2, hlm_B3),
       stars = c(0.01, 0.05, 0.10),
       include.ci = F,
       custom.header = list("Democratic Satisfaction" = 1:3,
                            "Democratic Evaluation" = 4:6),
       custom.coef.names = 
         c("Constant", "Performance: Narrowing Income Gaps", "Female (= 1)", "Age",
           "Age$^2$/100", "Education (4 = highest)", "Internet Usage (3 = highest)",
           "Performance: Managing the Economy", "Performance: Creating Jobs",
           "Performance: Keeping Prices Down", "Performance: Providing Health Services",
           "Performance: Addressing Educational Needs"),
       custom.note = "Entries are linear mixed-effects estimates with standard errors in parentheses.
       Missing data on individual predictors are imputed using MICE with default specifications in R (van Burren and Groothuis-Oudshoorn 2011).
       All significance tests are two-tailed with the following notations:
       $^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.1$.",
       caption = "Perceived Government Performance in Narrowing Income Gaps, Democratic Satisfaction, and Democratic Evaluation in African and Middle Eastern Countries",
       fontsize = "small",
       label = "table:individual_cor_me1")

## Table 2: Panels C and D (mixed effects)
texreg(list(hlm_C1, hlm_C2, hlm_D1, hlm_D2),
       stars = c(0.01, 0.05, 0.10),
       include.ci = F,
       custom.header = list("Democratic Satisfaction" = 1:2,
                            "Democratic Evaluation" = 3:4),
       custom.coef.names = 
         c("Constant", "Perceived Fairness of Income Distribution", "Female (= 1)", "Age",
           "Age$^2$/100", "Education (4 = highest)", "Internet Usage (3 = highest)"),
       custom.note = "Entries are linear mixed-effects estimates with robust standard errors in parentheses.
       Covariates on perceived government performance are not controlled for as data are not available among Asian and Latin Amiercan respondents.
       Missing data on individual predictors are imputed using MICE with default specifications in R (van Burren and Groothuis-Oudshoorn 2011).
       All significance tests are two-tailed with the following notations:
       $^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.1$.",
       caption = "Perceived Fairness of Income Distribution, Democratic Satisfaction, and Democratic Evaluation in Asian and Latin American Countries",
       fontsize = "small",
       label = "table:individual_cor_me2")

### PCA of different performance variables ----
## Load the required packages
library(FactoMineR) # version 2.8
library(factoextra) # version 1.0.7

## Create a data frame (with performance variables only) for PCA
df_pca <- df %>%
  dplyr::filter(country_un == 12 | country_un == 72 | country_un == 108 |
                  country_un == 120 | country_un == 132 | country_un == 204 |
                  country_un == 266 | country_un == 275 | country_un == 288 |
                  country_un == 324 | country_un == 384 | country_un == 400 |
                  country_un == 404 | country_un == 422 | country_un == 426 |
                  country_un == 430 | country_un == 450 | country_un == 454 |
                  country_un == 466 | country_un == 480 | country_un == 504 |
                  country_un == 508 | country_un == 516 | country_un == 562 |
                  country_un == 566 | country_un == 678 | country_un == 686 |
                  country_un == 694 | country_un == 710 | country_un == 716 |
                  country_un == 748 | country_un == 768 | country_un == 788 |
                  country_un == 800 | country_un == 818 | country_un == 834 |
                  country_un == 854 | country_un == 894) %>%
  dplyr::select(perf_incgap, perf_econ, perf_job, perf_price, perf_health, perf_educ) %>% 
  rename(`Narrowing Income Gaps` = perf_incgap,
         `Managing the Economy` = perf_econ,
         `Creating Jobs` = perf_job,
         `Keeping Prices Down` = perf_price,
         `Providing Health Services` = perf_health,
         `Addressing Educational Needs` = perf_educ)

## Conduct the PCA and generate a scree plot (Figure S5)
data.pca <- prcomp(na.omit(df_pca))
fviz_eig(data.pca, addlabels = T,
         barfill = "grey", barcolor = "black",
         main = "", ylab = "Percentage of Explained Variances", ylim = c(0, 100),
         ggtheme = theme_bw())
ggsave("pca_scree.pdf", width = 6, height = 4)

## Generate loading plots (Figure S6)
PCA(na.omit(df_pca), ncp = 3, graph = T, axes = c(1, 3))
ggsave("pca_dim3a.pdf", width = 5, height = 5)
PCA(na.omit(df_pca), ncp = 3, graph = T, axes = c(2, 3))
ggsave("pca_dim3b.pdf", width = 5, height = 5)

### Within-country analysis (fixed effects) ----
## Democratic satisfaction and govt performance in narrowing income gaps
reg_A1 <- lm_robust(dem_satis ~ perf_incgap, 
                    data = df_hlm_A, 
                    fixed_effects = country_un)
reg_A2 <- lm_robust(dem_satis ~ perf_incgap + 
                      female + age + I(age^2/100) + educ + internet, 
                    data = df_hlm_A, 
                    fixed_effects = country_un)
reg_A3 <- lm_robust(dem_satis ~ perf_incgap + 
                      female + age + I(age^2/100) + educ + internet + 
                      perf_econ + perf_job + perf_price + perf_health + perf_educ, 
                    data = df_hlm_A, 
                    fixed_effects = country_un)

## Democratic evaluation and govt performance in narrowing income gaps
reg_B1 <- lm_robust(dem_eval ~ perf_incgap, 
                    data = df_hlm_B, 
                    fixed_effects = country_un)
reg_B2 <- lm_robust(dem_eval ~ perf_incgap + 
                      female + age + I(age^2/100) + educ + internet, 
                    data = df_hlm_B, 
                    fixed_effects = country_un)
reg_B3 <- lm_robust(dem_eval ~ perf_incgap + 
                      female + age + I(age^2/100) + educ + internet + 
                      perf_econ + perf_job + perf_price + perf_health + perf_educ, 
                    data = df_hlm_B, 
                    fixed_effects = country_un)

## Democratic satisfaction and perceived fairness of income distribution
reg_C1 <- lm_robust(dem_satis ~ fair_incdist, 
                    data = df_hlm_C, 
                    fixed_effects = country_un)
reg_C2 <- lm_robust(dem_satis ~ fair_incdist + 
                      female + age + I(age^2/100) + educ + internet, 
                    data = df_hlm_C,  
                    fixed_effects = country_un)

## Democratic evaluation and perceived fairness of income distribution
reg_D1 <- lm_robust(dem_eval ~ fair_incdist, 
                    data = df_hlm_D, 
                    fixed_effects = country_un)
reg_D2 <- lm_robust(dem_eval ~ fair_incdist + 
                      female + age + I(age^2/100) + educ + internet, 
                    data = df_hlm_D, 
                    fixed_effects = country_un)

## Table S1: Panels A and B (fixed effects)
texreg(list(reg_A1, reg_A2, reg_A3, reg_B1, reg_B2, reg_B3),
       stars = c(0.01, 0.05, 0.10),
       include.ci = F,
       custom.header = list("Democratic Satisfaction" = 1:3,
                            "Democratic Evaluation" = 4:6),
       custom.coef.names = 
         c("Performance: Narrowing Income Gaps", "Female (= 1)", "Age",
           "Age$^2$/100", "Education (4 = highest)", "Internet Usage (3 = highest)",
           "Performance: Managing the Economy", "Performance: Creating Jobs",
           "Performance: Keeping Prices Down", "Performance: Providing Health Services",
           "Performance: Addressing Educational Needs"),
       custom.note = "Entries are OLS estimates with robust standard errors in parentheses.
       Missing data on individual predictors are imputed using MICE with default specifications in R (van Burren and Groothuis-Oudshoorn 2011).
       Country fixed effects are not shown.
       All significance tests are two-tailed with the following notations:
       $^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.1$.",
       caption = "Perceived Government Performance in Narrowing Income Gaps, Democratic Satisfaction, and Democratic Evaluation in African and Middle Eastern Countries (Fixed Effects)",
       fontsize = "small",
       label = "table:individual_cor_fe1")

## Table S2: Panels C and D (fixed effects)
texreg(list(reg_C1, reg_C2, reg_D1, reg_D2),
       stars = c(0.01, 0.05, 0.10),
       include.ci = F,
       custom.header = list("Democratic Satisfaction" = 1:2,
                            "Democratic Evaluation" = 3:4),
       custom.coef.names = 
         c("Perceived Fairness of Income Distribution", "Female (= 1)", "Age",
           "Age$^2$/100", "Education (4 = highest)", "Internet Usage (3 = highest)"),
       custom.note = "Entries are OLS estimates with robust standard errors in parentheses.
       Covariates on perceived government performance are not controlled for as data are not available among Asian and Latin Amiercan respondents.
       Missing data on individual predictors are imputed using MICE with default specifications in R (van Burren and Groothuis-Oudshoorn 2011).
       Country fixed effects are not shown.
       All significance tests are two-tailed with the following notations:
       $^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.1$.",
       caption = "Perceived Fairness of Income Distribution, Democratic Satisfaction, and Democratic Evaluation in Asian and Latin American Countries (Fixed Effects)",
       fontsize = "small",
       label = "table:individual_cor_fe2")
